Topological methods for searching barriers and reaction paths 
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We present a family of algorithms for the fast determination of reaction paths and barriers in phase 
space and the computation of the corresponding rates. The method requires the reaction times be 
large compared to the microscopic time, irrespective of the origin — energetic, entropic, cooperative 
— of the timescale separation. It lends itself to temperature cycling as in simulated annealing and 
to activation-relaxation routines. The dynamics is ultimately based on supersymmetry methods 
used years ago to derive Morse theory. Thus, the formalism automatically incorporates all relevant 
topological information. 



A situation often met in chemical reactions, protein 
folding and nucleation transitions is that the evolution 
breaks into fast local relaxations and rare activation 
processes. This timescale separation can be a conse- 
quence of low temperature (intra- valley equilibration oc- 
curs in Tfast — 0(1), passage times are exponential 
Tactiv ~ exp(A_E/r)), or of cooperative dynamics: e.g. 
in a ferromagnet of size L the relaxation time within a 
state Tfast < (the time for the largest domain to col- 
lapse), and between phases TacUv ~ exp(ci'^~^). It can 
also be the result of of entropic barriers, flat but narrow 
energy canyons that the system is unlikely to transit. 

In all these knowledge of the actual connectiv- 

ity of paths and saddles gives a global picture of a sys- 
tem's behavior . Direct information on the barriers can 
also be used to model the activated processes, without 
the need to wait for them to occur spontaneously, thus 
allowing to simulate a reaction even faster than nature. 
Furthermore, activation/relaxation procedures 0, alter- 
nating ordinary relaxation with forced activation, offer an 
efficient way to sample phase-space. For these investiga- 
tions, several algorithms P, |3| have been implemented to 
locate saddle points in energy. In cases when the nature 
of barriers is qualitatively affected by entropy, one needs 
to study directly the transition probability — thus taking 
into account the multiplicity of paths through the ( "free 
energy") barrier. A well known method to do this is a 
MonteCarlo sampling of the trajectories with ends fixed 
in the initial and the target state, known as 'Transition 
Path Samphng' H, In practice, this means going 

from an 'open-ended' to a 'two-ended' procedure. 

The purpose of this paper is to present a family of 
methods for the determination of reaction distributions 
and barriers. It is based on a dynamics that converges 
to barriers and reaction paths, much in the same way 
that diffusion in a potential converges to energy minima 
and metastable states. This allows us to construct 'free 
energy barrier' algorithms that can however be imple- 
mented in an open-ended manner. 

The organization of this paper is as follows: in the first 
section we present the dynamics and we justify it in the 
second section. 



DIFFUSION DYNAMICS AND PATH SAMPLING 

Consider a Langevin process in an A^-dimensional en- 
ergy landscape E{xi, xm)'- 



^-E,+r],{t) 



(1) 



where ri(t) is a white noise of variance 2T. Here and in 
what follows Ei and Eij denote derivatives of E, di = 
and we adopt the summation convention. An equivalent 
way of stating is to consider the probability distribu- 
tion P{x, t), evolving according to: 



dP{x,t) 
dt 



~HppP{x,t) = -d,Mx) (2) 
-a,; {Td, + E{) (3) 



which defines the current Ji{x, t) = {Tdi + Ei) P{x, t). 

The problem is how to modify in order to devise a 
diffusive dynamics that targets saddle points and reaction 
paths. We shall show below that a set of N functions 
Ri{x) evolving with 



dRi{x, t) 
dt 



= -HppR,{x,t) - E,j{x)R.j{x,t) (4) 



yields, at times Tfast < t <^ TacUv a vector field Ri{t) — 
Ji{x) describing reaction currents between metastable 
states — which one depends on the initial conditions. 
In other words, ^ does for reaction paths what ^ does 
for states. 

Equation is not immediately suitable for practical 
computations. Just as the natural way to simulate ^ 
is to use a diffusion dynamics fQ), in order to have a 
practical approach we have to find the diffusion dynamics 
yielding 

Vector viralkers. We consider a population of inde- 
pendent walkers x°'{t) carrying a vector internal degree 
of freedom f"(t), an A^-dimensional normalized vector. 
The dynamics consists of three types of transformations. 
Putting A{x°-) = —vfEijVj, they are : 

• Ordinary diffusion of the positions x°'{t) as in 

• Length preserving rotation of the vector with a po- 
sition-dependent rate: ^ = —Eijv'- {t) — A{x°-)vf ] 
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FIG. 1: Vector walkers describing a reaction path. The po- 
tential and the temperature are the ones ofRef. 0. There 
are two minima at (±\/5/4, 0), two saddles at (0, ±1) and a 
maximum at (0, 0). 

• Cloning with rate A{x'^) (if A{x°-) > 0), or death 
with rate —A{x'^) (if A{x°') < 0). Clones are born 
in the same place with the same vector jlH . 

The rotation-cloning steps play the role of singling out 
the least stable direction, somewhat like in eigenvector- 
following and related methods 0, Q • After a transient 
Tfast < t <^ Tactiv , the systcm tends to a regime in which 
particles organize along trails: the reaction paths (Fig. 

The distribution is almost stationary, there are re- 
gions of birth and emigration (near saddles) and regions 
of immigration and death (within the states). There is 
only a small global death rate of the order of the inverse 
reaction time along the specific path. 

Denoting F{x, v, t) the density of walkers, the current 
distribution is calculated from the average Ri{x,t) = 
J dv V Fix, V, t). One can check that the Ri then satisfy 
(@J), as desired. This algorithm can be seen as a Diffusion 
Monte Carlo as applied to ^ (see Refs. |1| and 0). 

Once an unnormalizcd transition current J is known, 
the reaction time is obtained for example through: 

_ fdx enj\' 

J dxePE{dWjy 

where /3 = ^. The barrier resistance e^^\J\'^ is large 
in (and serves to indicate) barriers, while e'^^(divJ)^ is 
concentrated and essentially constant within the states 
involved in the passage P|. In Fig. |21we show the barrier 
resistance, as obtained from the present algorithm, for a 
purely entropic barrier. 

Let us conclude by mentioning that sign problems are 
not a priori excluded here. They consist in the appear- 
ance of neighboring walkers with opposite vectors. For 
very low temperatures, and in particular for the evalu- 
ation of saddle-points, one can show 8] that this poses 



no problem. Even at finite temperatures walker cancella- 
tions turn out not to be serious. However, one can always 
decimate the walkers, lumping them together according 
to the mean local orientation. 

Transition path sampling with saddle-seeking 
paths. One can also use these ideas to improve tran- 
sition path sampling |^ Q . This family of techniques 
involves sampling over trajectories, weighted with the 
appropriate action. In order to obtain information on 
reaction paths, it is necessary to force the ends of the 
trajectory to be at either side of a barrier — otherwise 
the trajectory collapses to a globular shape around one 
minimum. 

If we wish the trajectory to seek for saddles, we need 
to add the effect of the second term on the rhs of Q . A 
simple calculation shows that this involves adding to the 
usual action a Lyapunov term equal to the log arithm of 
the largest eigenvalue of U(t) defined from [1^: 

Uij = —EikUkj (6) 

(Adding in general the fc-th eigenvalue makes the path 
seek the saddles with k unstable directions). 

One can show that the path obtained with the Lya- 
punov term is concentrated with a density proportional 
to the barrier resistance e'^-^jjp (cfr. Eq. jSJ): thanks 
to this modification, if the trajectory is not constrained 
in any way, it will spend its time near a barrier (even 
if entropic), and not in a state (see Fig. I^J, although 
the value of the action remains the same as for a pure 
Langevin weight. Thus the sampling can take better ad- 
vantage of the regions that are difficult for the passage, 
also eliminating the problem of the large fluctuations as- 
sociated with the jump time. If instead one end is tied to 
a minimum, the trajectory leaves the minimum and col- 
lapses to a globular shape around a nearby saddle (see 
Fig. 13): this can be used as a strategy to explore exit 
paths from a known configuration. The transition times 
in this setting can be calculated adapting the methods in 

Is]. 

THEORETICAL BACKGROUND 



3 




FIG. 3: Above: two-ended transition path sampling. Left: 
Langevin weight, right: Langevin + Lyapunov weight. Below: 
Langevin + Lyapunov action with only one end fixed while 
the other end is around a saddle. The potential is as in Fig. 

mi 



To show these results, and to place them in a wider 
context, we introduce the supersymmetric formaUsm. 
Let us start by artificially 'completing the square' in © 
as follows: we introduce the N fermion creation and an- 
nihilation operators Oi and aj , with anticommutation re- 



lations [ai,aj]-|- = Sij, and define charges as: 



Q 



(7) 



Both operators Q and Q commute with H, are nilpotent 
g2 ^ Q2 ^ 0, and: 

1 



H = -{Q + QY = H, 



(8) 



We consider a wave function ji/j) in the enlarged space 
and: 



dt 



(9) 



which yields a separate evolution for each fermion num- 
ber subspace. The zero-fermion subspace is just the 
original space of probability distributions and © is the 
Fokker-Planck equation © for them. The one-fermion 
subspace is made of iV-component functions Ri{x)al\) 
(I) the fermion vacuum) and © reduces to Q for them. 
More generally, the /c-fermion subspace is spanned by 
N\/k\{N -k)\ functions. 

As is well known, iJpp can be taken to a Hermitian 
form iJ^p = e''^/2^ppe-'''^/2, and this is also true of H: 



1 



T 



The eigenvalue equations are: 



-Ej, 



EijO^-ai 



(10) 
(11) 
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FIG. 4: Morse theory. The arrows indicate the action of Q. 
The gap in the spectrum reflects timescale separation. The 
numbers between brackets indicate the number of states below 
the gap for each fermion number. The Morse inequalities can 
be read from the picture. 



where 



-I3E/2 



(12) 



Now, it is easy to see that applying Q to any eigenvec- 
tor \'tl}^) we get either a degenerate eigenvector with one 
less fermion, or zero. Similarly, applying Q we get either 
a degenerate eigenvector with one more fermion or zero. 
Clearly, both Q and Q annihilate the Gibbs measure. 
The spectrum is thus organized as in Fig 21 

The subspaces with non-zero fermion numbers are at 
this stage artificial. The main result we need here is that 
small- eigenvalue (right) eigenvectors with one fermion 
encode the transition currents between states. On the 
other hand, small- eigenvalue k-fermion eigenvectors in 
the Hermitian base are peaked in all the saddle points 
with exactly k unstable directions. 

States. Before discussing transition distributions, 
we need to specify what we understand in general by 
'states'. In cases in which there is a timescale separation 
Tfast ^ Tactiv, whatcvcr its origin, the Fokker-Planck 
spectrum has a gap ~ r 



fast 



.(j^. Suppose there are if 
eigenvalues 'below the gap'. An important and intuitive 
result due to Gaveau and Schulman |9j is that using linear 
combinations of them we can construct K distributions 
that are essentially Gibbsean in K disjoint regions, and 
elsewhere negligible. In other words, metastable states 
(up to a lifetime TacUv) are linear combinations of states 
below the gap, and conversely. All these statements be- 
come sharper the larger the gap. The case of small tem- 
peratures is emblematic: the 'states' are simply narrow 
distributions (typically of width ^/T) peaked around lo- 
cal minima. This can be shown immediately from the 
low T expansion of (|10|l . 

Reaction Paths. Let us now sketch a proof of the 
fact that one-fermion eigenvectors below the gap give the 
reaction-path distribution between states, and their times. 
This means that Eq. Q on times larger than Tfast gives 
the reaction paths, and this in times much smaller than 
the passage-times Tactiv themselves. 
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Consider for simplicity the low-temperature situation, 
in which there are K minima separated by energy bar- 
riers, and a Langevin evolution starting from each mini- 
mum with a certain probability. If the time is of the or- 
der of the inverse of the largest passage time between any 
of the K states, there will be a current leaking through 
this saddle, and negligible flowing through all other paths 
(which take exponentially longer times). Next, consider 
times enough to allow the second fastest passage, there 
will be current flowing only through it: the faster re- 
action has already taken place (the states involved have 
mutually thermalized) and the slowest reactions-currents 
are exponentially smaller. If we continue this reason- 
ing K times, each time allowing for a new activation, 
we obtain a sequence of K different linear combinations 
(j)K{x), (I)i{x) of the states below the gap, the last one 
(j)i{x) the Gibbs measure. 

Consider next the K —1 onc-fcrmion right eigenvectors 
obtained acting with Q on the (p-zix), (j>K{x): 

Q<j)^ = (m + i?,Oal0^ = Jn=^)4\) (13) 

From the preceding argument, they represent exactly 
K — 1 independent passage current distributions. Hence, 
in general each individual one-fermion eigenvector below 
the gap obtained as in represents a passage current 
distribution with sources and sinks in the states (them- 
selves defined as above). In the low temperature case, 
the distribution is non-negligible along a tube following 
ascending and descending gradient lines, passing through 
a barrier. The direction of Ji{x) is parallel to the tube. 
The converse is also true: each passage between any two 
states is a linear superposition of these K — 1 currents. 

In Fig. 21 we see that there are also one-fermion eigen- 
vectors that are not the result of acting with Q as in 
(|13|l . They are characterized by being annihilated by 
Q, and this implies that they correspond to divergence- 
less current distributions. They are 'tours', for example 
an activated process leading from a minimum to itself 
through a saddle (think for example of a tilted Mexi- 
can hat). They are automatically distinguished in this 
formalism: for example (jSJ yields infinite timescales for 
them. 

Hermitian base: saddles and Morse Theory. 

One way to convince oneself that the formalism is a nat- 
ural one is to see how it yields quite simply the topo- 
logical relations between saddles (Morse theory). Let us 
review briefly the derivation in [lOj, specializing to the 
case of the topology of ordinary flat space and bounded 
systems. Consider the low temperature (harmonic) spec- 
trum of (|10|l : it is easy to show that fc-fcrmion eigen- 



vectors that have zero eigenvalue to leading order, are 
narrow Gaussians sitting on saddles with k unstable di- 
rections. This one-to-one relation between eigenfunctions 
'below the gap' in the Hermitian base and saddles implies, 
as can be seen from Fig. 0|that there are relations be- 
tween the numbers of saddles with successive numbers of 
unstable directions: (1 + ni), (ni -1-77,2), ... , ("-w-i -t-^Tv). 
The fact that Ua > constitute the strong Morse inequal- 
ities for our case. 

CONCLUSIONS 

We described a family of methods to find transition 
paths and saddle points that, to the best of our knowl- 
edge, are quite different from the standards ones. The 
statement of the algorithms is simple, but the inspira- 
tion for them - and indeed the proof that they work - 
is based on the supersymmetric quantum mechanics con- 
struction, and its connection to Morse theory. Although 
we could not provide here many details of the derivations, 
and we have not yet performed tests in really hard cases, 
we hope we have argued convincingly that the tools are 
quite relevant for the problem 0. 
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